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A Finite Element-Boundary Integral Method 
for Cavities in a Circular Cylinder 

Leo C. Kempel and John L. Volakis 
December 22, 1992 


Abstract 

Conformal antenna arrays offer many cost and weight advantages 
over conventional antenna systems. However, due to a lack of rigorous 
mathematical models for conformal antenna arrays, antenna design- 
ers result to measurements and planar antenna concepts for designing 
non-planar conformal antennas. Recently, we have found the finite 
element-boundary integral method to be very successful in modeling 
large planar arrays of arbitrary composition in a metallic plane. Here- 
with, we shall extend this formulation to conformal arrays on large 
metallic cylinders. In this report, we develop the mathematical for- 
mulation. In particular we discuss the shape functions, the resulting 
finite elements and the boundary integral equations, and the solution 
of the conformal finite element-boundary integral system. Some valid- 
ation results are presented and we further show how this formulation 
can be applied with minimal computational and memory resources. 
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1 Introduction 


Conformal antenna arrays are attractive for aircraft, spacecraft, and land 
vehicle applications since these systems possess low weight, flexibility, and 
cost advantages over conventional antennas. The majority of previous de- 
velopments in non-planar conformal antennas has been conducted experi- 
mentally due to a lack of rigorous analysis techniques. Various approximate 
analysis techniques are restricted in many respects, including accuracy and 
element shape, and are based on planar antenna models. 

Recently, we have found that the finite element-boundary integral (FEM- 
BI) method can be successfully employed for the analysis of large planar 
arrays of arbitrary composition [1]. The resulting system is sparse due to the 
local nature of the finite element method whereas the boundary integral is 
convolutional, thus ensuring an O(N) memory demand for the entire system. 

In this report we will extend the FEM-BI formulation for aperture an- 
tennas conformal to a cylindrical metallic surface. Both the radiation and 
scattering problems will be developed in the context of the FEM-BI method. 
In contrast to the planar aperture array, the implementation of the cylindric- 
al^ conformal array requires shell-shaped elements rather than bricks, and 
the required external Green’s function is that of the circular perfectly con- 
ducting cylinder. In its exact form this Green’s function is an infinite series 
which must be evaluated efficiently and must also be put in a convenient 
convolutional form for storage minimization. 

This report presents the FEM-BI formulation, the appropriate cylindrical 
shell elements, and the system evaluation strategy which will maintain low 
memory and computational load. The cylindrical elements will be chosen 
divergenceless while maintaining excellent geometrical fidelity. These ele- 
ments will be analogous to the bricks used by Jin and Volakis [1]. Efficient 
asymptotic evaluations of the cylindrical Green’s function will be discussed. 
In addition, the numerical implementation of the FEM-BI method will be 
presented along with some numerical results for validation purposes. These 
results will illustrate the accuracy of this formulation for the calculation of 
cavity eigenvalues and Radar Cross Section (RCS) calculations. Antenna 
pattern calculations will be presented in a future report. 
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2 FEM-BI Formulation 

Consider the configuration illustrated in figure 1 where a cavity is recessed 
in an infinite, circular metallic cylinder. The cavity walls are assumed to 
coincide with either constant p-,4>- or z-planes and we allow the possibility of 
radiating elements on the surface of or within the inhomogeneous substrate 
which fills the cavity. 

The FEM-BI formulation [1] allows the determination of the electric field 
present in the cavity due to interior or exterior sources. This method utilizes 
the finite element method to formulate the interior field and yield a sparse 
system of equations. This system is coupled to a second set of equations 
generated by enforcing a boundary integral equation on the cavity aperture. 
Although in general the boundary integral system of equations if full, a ju- 
dicious choice of boundary elements allows 0(N) storage to be maintained. 
The scattered or radiated field is readily calculated from the aperture fields. 

The development of the FEM-BI system begins with the vector wave 
equation valid in the interior of the cavity where we also allow the possibility 
of interior electric (J tni ) and/or magnetic ( M tnt ) sources. Specifically, we 
have 


V x 


V x Ejp^z) 
p r (p, <f>,z) 


k Ur{p, <f>,z)E(p, <f>,z) = 


- jk 0 Z 0 J ,nt (p,(f>,z ) + V x 


Af <n W,z) 

p r (p, 4>,z ) 


(i) 


where E(p,<f>,z ) is the total electric field, e r (/?, </>, z) and p r {p, <}>, z) is the 
relative permittivity and permeability of the substrate, k Q is the free-space 
wavenumber and Z 0 is the free-space intrinsic impedance. Unless otherwise 
noted, the time dependency is assumed and suppressed. To generate a 

system of equations from (1) we apply the method of weighted residuals 

/, { v * 

-kle r {p,<t>,z)E(p,<f>,z) ■ Wi(p,4>,z) \pdpd<f>dz = 


V x Ejp^z) 
p T {p,<t>,z) 


Wi{p, <f>, z) 
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/„ 


J U Y 

M int (p,4>,z) 

\ V x 

Pr{p,4>,z) 


jk 0 Z 0 J' n (p, <f>, z) > • W t (p, <f>, z)pdpd<pdz 


( 2 ) 


where W{(p, z) are subdomain vector weight functions to be specified. The 
substrate is discretized using cylindrical-shell volume elements such as the one 
shown in figure 2. The vector weight function coefficients represent the field 
at the edges of these elements and the integration volume ( V t ) corresponds 
to the elements which possess the i th edge. The forcing function due to the 
interior impressed sources (J tnt ,M ,nt ) is given by 




fl/ int (p, <f>, z) 
p r (p,(f>,z) 


jfc 0 Z 0 J (p, <j>, 


Wi(p, <f>, z)p dp dq b dz 


( 3 ) 


We recognize (2) as the weak form of the wave equation and upon ap- 
plication of a standard vector identity and the divergence theorem [2] (2) 
becomes 

f f V x E(p, (f>, z) • V x Wj(p, <f>, z) 

JVi\ fi r (p,<f>,z) 

-kle r (p,<f>,z)£(p,(j>,z) ■ Wi(p , <j>,z)^pdpd<j>dz 

-jk 0 Z 0 <f h{p , <f>, z) x H(p, <f>, z ) • Wi(p, (f>, z)dS = f\ nt (4) 

JSi 

with h(p,<f>,z) indicating the outward directed unit normal for the surfaces 
of the elements associated with the i th edge, 5, is the total surface of those 
elements, and H(p, <f>, z) is the total magnetic field. It can be shown that 
the surface integral in (4) vanishes for all surfaces except where the surface 
corresponds to the aperture. Thus, S, is the subdivision of the aperture 
surface associated with the elements which possess the i th edge. 

Unfortunately, (4) requires the determination of H(a,<f>,z) over the aper- 
ture which increases the required number of unknowns. To eliminate H(a, 4>,z) 
from (4) we use the integral representation of the magnetic field due to surface 
fields 

H(a,<f>,z) — H'(a, 4>, z) + H r (a, <f>, z) + 
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j¥ 0 k 0 a p(a , z) x £(a, z') • (7 0 (a, <£, dz 

+a J p{a , 2*) x H(a, <f>, z) • V x 6*0(0, <}>, z)d<j> dz (5) 

with 4> = 4>-<t>',z = z-z and 5 encompassing the entire surface of the 
cylinder. The first two terms are the incident field and the scattered field 
from the cylinder with the aperture removed. We shall denote these two 
terms by the symbol H cyl (a,<f),z). 

The integral equation (5) cannot be readily coupled with (4) because 
it involves both the tangential electric and magnetic fields at the aperture. 
However, we are free to choose a dyadic Green’s function which satisfies the 
radiation condition and the additional Neumann surface boundary condition 

VxC n2 (M,f) = 0. (6) 

With this condition, the second integral in (5) vanishes and the integration 
surface need only extend over the aperture surface. The dyadic Green’s 
function which satisfies (6) is usually termed the magnetic dyadic Green’s 
function of the second kind [3]. We may now write (5) as 

H(a,^z) = H cyl (a,<j>,z) 

+jY 0 k 0 a I p(a, <(>, z) x E(a, <f> , z) • G m2 (a, <f>, z)d<p dz (7) 
J s a 

where S a indicates integration over the entire aperture. 

Upon inserting (7) into (4), we obtain 

j f V x E(p, MWx W,(/>, <j>, z) 

-lcle T (p,<t>,z)E(p,<l>,z) ■ Wi(p,(f>,z)^pdpd<t>dz 

-\-(k 0 a) 2 / \Wi{a,<i>, z) ■ p(a,<f>,z) x G m2 {a,<t>,z) x p(a,<f> , z')- 

J ^> t J s 

E(a , <j>\z')\ d<j>' dz'dfdz = fl ni + ft xt ( 8 ) 

in which a standard dyadic identity has been used and the forcing function 
due to exterior sources is 

fi Xt = jZ 0 k 0 a j IT, (a, <f> , z') • p(a, <f> , z) x H cyl (a,<f>\z')d<f>' dz (9) 
Js, 
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We recognize that (8) has unknowns corresponding to the electric field within 
the cavity and on the aperture. 

Following the principles of Galerkin’s method, the electric field in (8) is 
expanded in terms of the same vector subdomain weight functions as was 
used for testing, i.e., 

= ZE,W,{p,i,z). (10) 

J=1 

In this expansion (interior -fi aperture edges) A 1 is the total number of un- 
knowns or edge fields. Since the vector wave equation (1) requires the electric 
field to be divergenceless, the weighting functions used in (10) must also be 
divergenceless to avoid the use of a penalty term [1] which is required if 
this condition is not met (e.g. node-based elements [4]). By necessity the 
aperture field takes the form 

E{a,<j>,z) = ^E 3 6 a {j)W } (a,<}>,z) (11) 

j'=i 


where 


6 a (j) = 1 if j € aperture edge 

= 0 else (12) 

and the vector weighting function is the same as that used in (10) when 
evaluated on the surface. Combining (8), (10) and (11) we obtain the FEM- 
BI system 

f f V x !Fj(p, <^>, z) • V x f-F,(/), <ft, z) 

Jv,\ p r (p,<f),z) 

- klt T {p,<f>,z)Wj(p,<f>,z ) • W,(p,4>,z)^pdpd4>dz 

+(k 0 a) 2 6 a (j)8 a (i) [ f [W i (a,^z)-p(a^,z)x 
J S t J s 

V m2 (a, l z) x p(a , 4>\ z) • Wj{p, <t>, z)} d<f>' dz'dcf>dz = // nt + fi* (13) 

in which Si and Sj represent the surface elements bordered by the i th and j th 
edges, respectively. Below we discuss the specific vector weighting functions 
for the shell element shown in figure 2. 
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3 Vector Weight Functions 

Traditional finite elements associate the unknown field coefficients (Ej) with 
the node field. However, it has been found that such elements enforce non- 
physical field continuity at element boundaries since both the tangential and 
normal field components are continuous which leads to spurious modes that 
must be suppressed with a penalty function [4]. However, edge-based ele- 
ments have only tangential continuity and are therefore better suited for 
electromagnetics applications where jump discontinuities in the normal com- 
ponent is permitted. In addition, edge-based functions avoid an explicit spe- 
cification of the fields at corners where the field may be singular [5]. Various 
volume element shapes have been used for finite element applications. Two 
of the most popular are: bricks [5] and tetrahedra [6, 7]. Although the later 
element is quite versatile, we shall develop a cylindrical shell edge-based ele- 
ments since we require geometrical fidelity which is not available with bricks 
or tetrahedra. 

Bricks are the natural element for discretization of a volume which is 
defined by portions of the Cartesian planes. Likewise, cylindrical shell ele- 
ments are superbly suited for cavities where the six faces are defined by 
portions of the constant /?-, <£-, and z-planes. Thus, the cylindrical shall ele- 
ments will be analogous to brick elements used by Jin and Volakis [5]. By 
inspection, the shell element in figure 2 is comprised of twelve edges which 
will constitute the vector weight functions. These weight functions should 
satisfy the following properties: 

1. subdomain (finite element) 

2. || Wj(p, <f>, z) ||= 1 if ( p,<f>,z ) € j th edge 

3. || Wj{p,<f>,z) ||= 0 if ( p,<t>,z ) <E any edge || j th edge 

4. V-W j (p,<f>,z) = 0 1 


l lVj(p y <f> y 2 ) will only satisfy this requirement within the volume of the element. These 
weighting functions introduce artificial charges on the faces of the element and are not 
divergenceless at element interfaces. This is required since these elements do not guarantee 
normal field continuity across the element faces. 
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These conditions are met by the following set of weight vectors 


p *4 

W l 2 \P, <t>, 2 ) = Wp{P, <Pr, 2 t , + ), 

*) W p{^p , (j), Z (f> r , Zb, ), 

^(p, <£, 2) = W^(/>, <f>, Z\ Pb, •, 2 t , +), 

WWp, </>, 2) = $*(p, <£, z; pt, z b , -), 


W^15(/>, 4>, 2 ) = W*(/>, <£, 2? P6, ^r, •, + ), 

W^sCp, <£, 2:) = W g (p, <j>, z; Pb, <f>t, 


W43 (p, <f>, z ) = W p (p, 4 > , 2; •, <£/, 2 t , -) 

U r 87 (/> ? <t>, - ) = IT P (/), <£, z\ *, (f>i , z t , + ) 

fT 23 (^, <£, *) = <£, z; />a, sZf, -) 

M / 67(/^i W Z » Pa 5 ’ 1 ^ 6 ? T) 

U 26 (p? ^ * ) h s(/^i /^oi ’? ) 

II 3?(/h 2) IE z (/9, 4*1 Z) Pa i *fih ’ ? T ) 


(14) 


where IE/;. denote the edge which is defined by local nodes (l,k) which is shown 
in figure 2. Each vector weight is represented by one of three fundamental 
vectors (one for each orthogonal direction) 


W p (p,4>,z\p, lz,s) 

~ > P' $*> ~ f ^ ) 


SPa (<j>-4>){z ~ Z) 
ah n 1 


-{p - p)(z - z)<fi 


W z (p,<t>,z',P,4>,z,s) 

with i = p b — p a , a = (fi r — (fit, and h = z t — z b . The element parameters 
(Pa,Pb> <fih <fir, Zb, z<) are shown in figure 2. We observe that all four of the 
requirements listed above are met by (14) with the understanding that these 
weight vectors are non-zero only within the element. Interestingly, in order 
to satisfy the divergenceless requirement by forcing each component of the 
element field to be divergenceless, we need the - term in the ^-directed 
weight (15). Of course, for very large radius cylinders and small elements, 
the curvature of these cylindrical shell elements decreases resulting in weight 
vectors which are functionally similar to the bricks used by Jin and Volakis 
[5], We shall now specify the appropriate dyadic Green’s function which will 
allow us, in conjunction with the vector weight functions given above, to fully 
specify the FEM-BI system (13). 


ta 


(p - p){<t> - 4>) l * 


(15) 
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4 Dyadic Green’s Function 

The boundary integral present in (13) requires numerical computation and 
an efficient evaluation of the dyadic Green’s function, G m ?(a,4>, 2), is crucial. 
This dyadic Green’s function may be derived in its eigenfunction form from 
a set of cylindrical vector wave functions using the procedure espoused by 
Tai [8]. The resulting expression contains an infinite sum of angular eigen- 
values and an infinite integration over the axial eigenvalues. It is well known 
that the angular eigenvalue expansion for cylinders exhibits extremely poor 
convergence characteristics for large radius cylinders [9]. We must therefore 
resort to approximate evaluations of the Green’s function which do not pos- 
sess the prohibitive computational demand of the exact modal expressions. 

As described above, one may readily derive a dyadic Green’s function 
which satisfies both the radiation condition (by virtue of the chosen wave 
functions) and the Neumann boundary condition on the surface of the cylin- 
der (6). The resulting modal expression is 

G m2 (a, <f>, z) — d»^G^(a, </!>, z) -f [<£r -f z<}> ] G 4>z (a, <j>, z) + zzG zz {a , 4>, z ) 

(16) 


where unprimed unit vectors are functions of the observation point (a,4>,z) 
while primed unit vectors are functions of the source point ,z ) and 


G zz (aJ,z) = 


l y' f°° ( 1 j(n»-fc«g)^r 

00 UJ 7 H '<%) 


± =\ _ 1 y' f°° ( n ^ z \ i(n<t>-k z z) 

G («.*,) - wJL'U Via-,) HP l( 7) e 


- 1 _ cc _ roc ] 


K m h) 


_ ( ok. V fff'l 

\hakj h'P 


n?Hj) 

(7) 


!) dk 


( 17 ) 


with k p — yjkl — k%, 7 = k p a, a is the radius of the cylinder, n denotes the 
angular mode number, and k z is the axial eigenvalue. As stated previously, 
these expressions exact a high computational demand especially when the 
source and observation point coalesce. We therefore want to recast (17) 
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in some other approximate form which improve as the radius of the cylinder 
increases and the distance between the source and observation points becomes 
large. 

Several research groups have developed in the past such evaluations of 
(17) and the most successful are attributed to Pathak [10], Boersma and Lee 
[11], and Bird [12, 13]. All three approximations utilize Watson’s transform- 
ation to convert the slowly converging angular eigenvalue representation (17) 
into a rapidly converging series over the radial eigenvalue. The Hankel func- 
tion and its derivative are approximated using either the Debye expansions 
or the uniform asymptotic expansion [14] and finally the axial eigenvalue 
integral is evaluated via the steepest decent method. These three formulas 
differ only in the level of approximation offered. The three formula are listed 
in tabular form by Bird [12] where he compared the accuracy of each. In 
developing these formula, the following quantities are necessary in describing 
the on-surface ray contributions: 

1. Geodesic path length fs = \J (a<I>) + 


2. Geodesic trajectory [0 = tan 1 ) 


3. 


Approximation order 



k 0 s 


) 


and in these ^ or = 2?r — </> depending on which path as illustrated in 
figure 3 (the short or long) is used. 

Boersma and Lee [11] employed approximate expressions from a Debye 
approximation of the Hankel functions and a first-order evaluation of the 
axial wavenumber integral to obtain 


G zz (a, <f>, z ) ~ “^ e ~ jfcoS { ( cos2 ° + 9(1 “ ?)(2 - 3 cos 2 e)) v(/3) 

-qfiv\(3) |j!co.s 2 0 + q (^- + | ian 2 6 - -^cos 2 tf) 

G* 2 (a, <fi,z) ~ 7 — : -qe~^ k ° 5 sin 0 cos 0< (1 — 3g(l — q)) v((3) 

2tt ( 

+<10v'(l3) ^tan 2 0 - 1 + q (~- - jLec 2 #) | 
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G <t><t> {a, <f>, z) ~ jfc ° S j ( stn2<? + ?(1 - ?)(2 - 3sm 2 0)) u(£) 

+<? [sec 2 0 («(/?) - «(/?)) + (3 (i>\/3) fyan 2 0 - ^*in 2 0 

u'{(3)^sec 2 e + v'(0)q (J - 2- tan 2 9 + ^p«n 2 0^ | 

(18) 

Pathak’s expressions are based on a uniform asymptotic expansions of the 
Hankel functions and a first-order evaluation of the axial wavenumber integ- 
ral. His results are 


G n M,z) 

G+ Z {a,$,z) 

G**(a,lz) 


— (cos 2 0 + q( 1 - q)( 2 — 3cos 2 #)) n(/3)j 
^^-qe~ jkoS sin 6 cos (1 — 3q(l — <7))t’(/?)j 

— ^-qe~* koS \^ (sin 2 0 + < 7(1 — q){2 — 3sin 2 0)j t;(/3) 
+<7 [sefO (u(0) - v(0)j\ | 


(19) 


Finally, Bird [12] used higher order terms from the uniform asymptotic for- 
mulas and a second-order evaluation of the integral to obtain 


G zz (aJ,z) 


G* z (a,lz) 


~^~Q e • 7fc ° 3 | ( cos 2 6 + ^(l — q)(2 — 3 cos 2 0)') v(/3 ) 

+q (—sin 2 0 — — 'j v(0) + ( —-sin 2 o\ V\([3) 

v L V72 247 V60 36 J y ’ 

{^sm 2 0 + -t) v 2 (I3) + J -f3v { 0 1] (/3) J 

-qe~i koS sin 0 cos o\ (1 — 3^r( 1 —q))v(/3) 

2 71 


+9 


IstcH - H'l 


( 17 90 \ 

- - —sec 2 0j Vl (/3) 


12 J 
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(■ J^ sec2 ° - V M + ^/3sec 2 0v { o l) (/3) | 

z) ~ -^~qe~ jkol> ^sin 2 0 + q(l — q)(2 — 3sin 2 6)^ v(P) 

+q {u(/3) — v{fl))sec 2 0 + (^-tan 2 6 — — sm 2 $) v(fl) 

\ 1 7 / 4-( / 

+ _ S ian2 ^) vi ^ 

+ (l5 tan2 ° ~ ^l sin2 °) l ’ 2 ^ + i ) P tan7&v ’o X) (^) | (20) 

2 

where /? = A:s [ ^p ~] 3 . In (18)-(20), u(/3) and v™(/3) represent the soft and 
hard surface Fock functions, respectively, while u(/3) and v'(fl) denote the 
first derivative of the zeroth order surface Fock functions. These functions are 
characteristic of the creeping waves on a circular cylinder and are discussed 
in detail by Logan [15]. 

We will consider only the two direct path contributions from the creeping 
wave series which are shown in figure 3. Accordingly, the ray parameters 
given above depend on which direct ray contribution is being computed. 
Bird [12] has concluded that his expressions (20), which include terms up to 
0{q~ 2 ) are the most accurate throughout the full range of ^ and z except 
when s is very small or the observation point is in the paraxial region of 
the source. In this region, the formula given by Boersma and Lee (18) are 
most accurate since the Debye approximation used for the Hankel function is 
well suited for the paraxial region and 0(q~ 3 ) terms have been retained. A 
notable exception to this conclusion occurs when the cylinder is sufficiently 
large while the aperture is small. In this situation, the aperture appears 
to be a hole in a planar metallic screen and the formula given by Pathak 
(19) turns out to be the most accurate. Inspection of (18)-(20) reveals all 
three approximate expressions include the free-space dyadic Green’s function 
terms up to 0(q~ 2 ) modulated by a surface Fock function whereas Bird’s 
as well as Boersma and Lee’s formula contain ’extra’ terms. These terms 
are not present in the case of a metallic screen and thus provide reduced 
accuracy when the aperture looks planar. However, for cavity configurations 
where curvature is considered, Bird’s conclusions should be valid. Therefore, 
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we shall use Boersma and Lee’s formulas (18) if |cos0| < 0.1 or s < 1A. 
Otherwise, Bird’s formulas (20) will be used. The derivation of (18), (19) and 
(20) is described in [11], [10] and [13], respectively. 


5 Numerical Implementation 


In the preceding sections of this report, we described the vector weight func- 
tions and dyadic Green’s function appropriate for the solution of (13). Below 
we now explicitly develop the admittance matrix and excitation vector entries 
necessary to solve (13) for the electric field. Once the unknown field is de- 
termined, engineering quantities such as cavity eigenvalues, RCS and antenna 
radiation patterns may be determined. 

We may write the N equations implied by (13) in terms of N unknowns(£^ ; ) 
as a matrix equation 



[G] [0] 
[ 0 ] [ 0 ] 



text 
J 3 

ft 11 


( 21 ) 


where [.4] is a very sparse N x N symmetric matrix, [ Q ] is an N ap x N ap 
square full matrix, the N ap x 1 column vector denotes field values on 

the aperture edges, the N v x 1 column vector { E ) nt } is field values on interior 
edges, and the total number of unknowns (N) is the sum of the free aperture 
edges (N ap ) and the free interior edges (N v ). A free edge is any edge which 
is not tangential to and on a metallic wall since such an edge is known to 
be associated with a null field due to the natural boundary condition on the 
metallic wall. Hence, it is not included as a system unknown. The N ap x 1 
column vector {/“*} is the forcing function attributed to an incident external 
field (9) or impressed sources on the aperture (3) while the iV„ x 1 column 
vector {/•"'} is due to the impressed sources within the substrate (3). 


5.1 Matrix Assembly 

The sparse nature of [4] is revealed upon inspection of (13). Since the edge 
vectors are non-zero only within the elements which possess that particular 
edge, the volume integral is identically zero unless both the source and testing 
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edges are within the same element. These matrix entries may be written as 


A - — 

71 u — 1 st 

dr 


k 2 f /^ ,J 

K o t r J st 


( 22 ) 


where we have assumed constant material properties within each element 
with e T and p r denoting the relative permittivity and permeability of the 
element. The subscripts (i,j) indicate the row and column of the matrix 
and are related to the global edge numbers. Since each edge may belong 
to more than one element, the matrix entry A tJ is computed by adding the 
contributions for each element which have both the i th and j th edges. The 
auxiliary terms 1 ^' J " and are defined by 


= J V x IF s (p, 4>, z\ fa, fa, Zj , Sj) ■ V x W t (p, <j>, z ; p { , fa, z t ,s,)p dp d<j> dz 
l£ )lJ = j W,(p,faz-,pj,4>j,Zj,Sj) -W t (p,<l),z-,pi,4>i,zi,s,)pdpd<i>dz (23) 


with (s,t) € { f) , <f>, z] indicating the direction of the source and test edges, re- 
spectively, while the weight vectors are given by (14). Since the fundamental 
weight vectors (15) have orthogonal orientations and (23) is symmetric with 
respect to the source and test vectors, 1^ need be determined for six com- 
binations of (s,t) and iffl need be determined for three such combinations. 
These integrals may be evaluated analytically upon performing the vector 
operations indicated in (23). 

In addition to [A], we must also evaluate the entries of [Q] which includes 
all the boundary integral contribution to the FEM-BI system (13). The 
entries of [Q] involve the integral 


G tJ = ( k 0 a ) 2 / / W t (a,<f>,z;p„<f>i,zi,si) ■ 

J 5 , JSj 

p(a, <p , z) x (7 m2 (a, <£, z) x p{a, <j>\ z) 
(a, <j>, z \ pj, <j>j, Zj,Sj)d<f>' dz'd<f)dz 


(24) 


which unfortunately cannot be evaluated analytically due to the presence 
of the Green’s function. From an inspection_of the approximate Green’s 
function formulas (18)-(20) we observe that (?m2( a , <f>, z) becomes singular 
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as the source and observation point coalesce. When these points^ arc close, 
a reasonable approach for large radius cylinders is to replace G m 2 (a,<^, z) 
with the metallic screen Green’s function whose singularity has been treated 
previously in the context of this integral [4]. Accordingly, we use the following 
boundary integral when the source and test edges are either in the same or 
some adjacent element 


Gjj = 2 (k 0 a) 2 j / 

JS t J Sj 

p{ a i ~) X G 0 {a,4>,z) x p(a,<f>',z') 
•W a (a, <f>, z'\ fa, <j>j, zj, s } )d<f>' dz' d(j)dz 


(25) 


where 


G 0 {a, <f>,z) = 


= VV 


e -jk 0 R 
4x R 


(26) 


with R = \r — r | indicating the distance between observation and source 
points and / = xx + yy + zz. Upon use of a common vector identity and the 
divergence theorem [4], we have 


Gij = 


{Kaf 


2tt 

,2 


f f m / t ii 

/ / W t (a,4>,z] Pi,4i,*i,5i) ■ 1 ■ W s (a,<f) ,z ] pj, (f>j, z^s^dj) dz d<j>dz 

Js t JSj 


Js Is V ' x W t (a,<f>,z\pi,<i>i,Zi,Si)} 


it- ii i e~ jkR i I 

V- [p(a,<f> ,z ) x W„{a,4> ,z'\pj^ h z h Sj)\ ——d<j)'dz'd<f>d. 


R 


(27) 


This form of the boundary integral may be readily evaluated even as R 
vanishes by employing the regularization procedure used by Jin and Volakis 
[4] . The asymptotic dyadic Green’s function used here consists of two direct 
interactions as shown in figure 3. We may use the planar approximation 
(26) for the short path contribution when the distance between the source 
and test elements is small. The asymptotic formulas are always used for 
the long path contributions, i.e. interactions between distant elements when 
curvature effects are likely to be prominent. 
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Explicit formulas for the computation of the admittance matrix entries 
have been given above which may be readily implemented. The [A] matrix 
may be considered the finite element portion of the system while \Q\ is the 
boundary integral contribution. The matrix [.4] is very sparse and its entries 
may be computed from analytical evaluations of the integrals (23). Although 
the [0] submatrix is generally full and its entries must be computed numeric- 
ally using for example two-point Gaussian quadrature, considerable storage 
and computation economy is possible by utilizing a CG-FFT type of solution 
procedure. The convolutional kernel of the dyadic Green’s function makes 
this possible (17). The computation of the forcing functions in (21) must 
now be addressed. 

5.2 Excitation 

Two principal types of excitation are found in the FEM-BI system (21). One 
is the excitation due to internal sources (possibly on the aperture) when the 
radiation is being considered. The entries for the resulting column vector 
{/?”*} is given by (3) and require the specification of the impressed currents 
( J ,nt , M' nt ). The second type of excitation is due to external sources which 
are non-zero for scattering problem computations. This excitation column 
vector {fj Xt } are given by (9) which requires the surface magnetic field ( H cyl ) 
when no aperture is present. This later excitation will now be examined 
assuming plane wave incidence. 

Circular cylinder scattering is one of the few exact solutions available in 
electromagnetics. Traditionally, the field scattered by a cylinder due to plane 
wave incidence is expressed in modal form. Assuming that the incident field 
is given by 

E' = Ve-M** 

IP = Y 0 {k l x 

= Y a \p‘ sin a cos ft - o' cos c. - i sin a sin ft] ft 1 " 1 """ 1 ’' •»(*-*>+««»*] 

(28) 

where (0,-,^,) indicate the direction of incidence, a is the polarization angle 
where e' = cosa0' + sina<^‘, and Y 0 is the free-space intrinsic admittance. 
The total surface field 

H cyl (a, <j>,z) = H i (a,<j>,z) + H 3 cyl (a,<t>,z) (29) 
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is given by 




P-Yc 


sm a 


7T (2 


jk 0 cos 8 t z 


^oo [Hn 2 \k 0 a sin 0,) J 


-2Y 0 


p jk o cos0 t z oo 


L 


7r k Q a sin 0, [Hi 2) (k 0 a sin 0 t ) 


cos a 


+ 


n 


sin a cos 0 % 


^ k 0 a sin 0 t Il} 2 \k 0 a sin 6{) J 


P jn(f +<£-<£.) 


( 30 ) 


These expressions may be computed by summing only a few terms of the 
series if k o as\n0, is small. However, as this parameter becomes large (e.g. 
for large a and 0, — > 90°), asymptotic evaluations similar to those used for 
the dyadic Green’s function must be applied to (30) to maintain computa- 
tional economy. Utilizing Watson’s transformation and Fock theory [9] in 
connection with (30) we find 


H c z yl ~ -Y 0 sin a sin 0p k ° cos6 ' z J2 e~ 3k ° a ™ e [sr {0) (m$ p )] * 

p=l 

Ilf ~ j2Y 0 cos a - — — r— r- e 3k ° cos8,z T e - 3k ° asin0 '*’> f/ (0) (m$ p )] * 
* k Q a sm6 t “f 1 J 

</ 0) (™*p) 

-J - /, -g (1) ( m ^p) 


p=l 


-Y 0 sin a cos 0 t e 3k ° cos e,z ^ ( — 1 ye ~ jk ° a sin 8 ' <1>p 

p=l 


m 


k Q a sin 


(31) 


in which <t, = f -(*-*), 4>, = (*-*)-f, m= [iaaiiai]’, and the 
indicates complex conjugation. The appropriate Fock functions are [15] 2 


/>(0 

/ (,) (0 


V5F 

jl 


L 

l 


r «h(0 

iO 


V tr>i(<) 


dt 


(32) 


2 Logan uses the e time dependency in the definition of these functions requiring 
the complex conjugation in (31) 
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where the Fock form of the Airy functions may be written in the more familiar 
Miller notation u-’i(f) = >/5r [#;(<) + jAi(t)] and the integration contour is 
given by Logan [15]. 

The asymptotic formulas are quite accurate compared to (30) except when 
<f> « <j > In this region, Goriainov [16] derived the following more accurate 
expressions from Fock theory 

Hf ~ 


Hf ~ 


with [15] 

G m ( 0 = 9 m (0^ 

F m (() = (34) 

These surface field expressions may now be used to efficiently calculate 
the entries of the column vector {/?*'} via a numerical evaluation of (9). In 
particular, the modal series (30) is used when k o as\n0i < 10 and either (31) 
or (33) for k 0 a sin 0, > 10 as appropriate. We note that the unit vectors in 
(9) and (28) are not the same since the former depend on the integration 
variable (<f>) while the latter depend on the incidence angle (<f>i). Thus, the 


-K sin a sin [/^(m#,)]’ 

+ e jk a as in* «.(*-*) [G(-m COS (<f> - <f>i))]‘ | 

J ko COS 0, 2 ( e -jkoa sin 6, j" y( (0) ( m ^ ) 


]2Y 0 cos a 


k 0 a sin Q t 


+c j* 0 ««n *«»<*-*,) [/r(_ mcoS ( < 6- &))]* | 
+Y 0 sin a cos ^. c i*oC0.* i *| c -jfc„«sin* i * I g ( 


m 


-J 


5 (1) (m$!) 


k 0 a sin 9 t 

jk 0 a sin 6 t cos (4>—<t > t ) 


G(— m cos ( <f> — <}>i)) 


-J 


m 


k Q a sin 0, 


—m cos (<f> — <f>i) 


■) 


(33) 
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(35) 


excitation integral (9) is 

f- xt - jZ 0 (ka) f cos (<f> - <f>i)Wi (a, <f> , z) • 

[i zHf(a , <f>\ z ) - 4> H c z yl (a, 4>\ z)\ d<f>' dz 

6 Engineering Quantities 

The FEM-BI system (21) has now been fully specified and either a direct or 
an iterative solution of (21) for the system unknowns {£,} may be applied. 
Once the expansion coefficients have been determined, we may calculate use- 
ful engineering quantities such as cavity eigenvalues, RCS and radiation pat- 
terns. 

6.1 Cavity Eigenvalues 

One of the uses of the Finite Element Method in electromagnetics which has 
become popular recently is the calculation of cavity resonances numerically. 
Often analytical determination of the cavity modes is not possible if the cavity 
has an odd shape or is filled with an inhomogeneous material. Chatterjee et 
«/. [7] have shown that edge-based tetrahedral finite elements are well-suited 
for computing cavity eigenvalues. These elements are able to model a wide 
range of cavity shapes and material filling. 

Likewise, the cylindrical shell elements introduced in this report are ideal 
for the calculation of modes present in any cylindrical cavity which does 
not intersect the z-axis. For eigenvalue calculations, the cavity aperture is 
covered with metal which closes the finite element mesh, implying [Q] = [0]. 
Therefore, the eigenvalues may be computed by setting the interior portion 
of the system (22) equal to zero and identifying the eigenvalues as A = k% 

— /jH (By) = A M? W ] {By) (36) 

.^r J 

We recognize (36) as a generalized eigenvalue problem where A denote the 
eigenvalues of the cavity. These eigenvalues may be determined directly or 
by exploiting the inherent matrix sparsity, Lanczos method may be used. 
Such a solution of (36) yield the resonant fields {Ej} (eigenvectors) and 
corresponding wavenumbers (eigenvalues). 
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6.2 Far-field Evaluation 


Two of the most important applications of the formulation presented in this 
report is the calculation of the cavity’s radar cross section (RCS) and the 
radiation pattern due to sources placed within or on the aperture of the 
cavity. In this section, we will describe far-field calculation due to equivalent 
sources on the aperture. This will entail a similar analysis as was performed 
for the excitation fields. 

We begin with the integral representation of the scattered magnetic field 
in terms of the aperture fields. We have, 

H s (r,9,<j> ) = jY 0 k 0 a J s G m2 {r,0,(f>]a,(f)',z') • 

[p(a, <j>\ z) x E(a, <f> , s')] d<f> dz (37) 


with (r ,0,<f>) indicating the observation point in spherical coordinates. When 
the observation point is very far from the cylinder, the dyadic Green’s func- 
tion has the asymptotic form 


b m 2( r i ? z ) 


c-jk 0 r 


k Q r 


G H 9(j) + G ez 6z + (38) 


where unprimed unit vectors are functions of the observation position while 
primed unit vectors are functions of the integration point in (37), The com- 
ponents of this asymptotic Green’s function are determined by a mode match- 
ing procedure 


G h 

G $z 

G M 


j 2k 0 cos 0 jkoCOS0z ‘ ^ 

■ e ^ TT '{ 


n 


(2?r) 2 (k 0 a sin 6) 2 “ H^\k 0 a sin 9) 

j 9 < 00 1 

J_ y ^ j k 0 cos Oz * ' 

(27t) 2 a 


gjnff +(<t>~<l> )) 


e 53 rr'(2) 


n= — 00 Hn [Jc o CL SIT! 0 ^) 


e i"(f +(<l>-<t> )) 


0 jk 0 cosSz 


1 


(27t) 2 a sin 9 


Hn\k 0 asm9) 


e i n (f +(4>-<t> )) 


(39) 


As one might expect, these series converge rather slowly for large k o asin0. 
They must therefore be recast in another form by employing Watson’s trans- 
formation and Fock theory as described before. We have, 


G h 


koC°s° e j ko cos ez' j2(-l)Pe-j k ° asin0 *r |g (0) (m3> p ) - j ff (1) (m$ p ) 

4 7r j L k Q a sin 0 
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— — — — — e^ k ° cos ^ z> Y" e~ jk ° a 
47r 


^ e -jk 0 a sin B* p L<0)( m $ 


2 , 2 

^ ^ ^jkpCosOz V"* jk 0 a sin (9<J> p 


2air sin 6 




where the appropriate Fock functions are given by (32). As expected in the 
region <f> fa <f>, the formulas attributed to Goriainov [16] 


G h ~ - - — 6 e jko cos6z ' i e~ jk ° a sin **» [<7 (0) (T7i^i)-i , m yi ff (1) (77?^ 1 )l 
4 7T [ k o as\u0 

eJ k 0 asm8costt-*') [ 6 -(0)(_ m CQS ^ j )) _ j ™ g 0)(_ m cos 

„ k 0 o> sin 0 


fc 0 Sin 6 ^jk 0 COS 0z' f e -jk c a sin 9*! + 


jk 0 a sin 8 cos 


<t> ) ^(°)(_ m cos — (j))) 


Q<t>4> ~ e jk 0 cos9z J e -j* o asin0<*i \ j \1 * 

2a7rsin0 1 >■ J 


g jfeoO sin 6 cos [ j p(»)(_ m cos _ <£'))]* 


are most useful. The far-zone scattered field can be computed numerically 
by using either (39), (40), or (41) in (37). For the scattering problem, RCS is 
most often the quantity of interest. Once the far-zone magnetic field (37) is 
computed, the RCS is given 


a(0,<f)) = lim 47rr 




7 Numerical Results 

The preceding section of this report presented several engineering quantities 
of interest: cavity eigenvalues (resonant modes), RCS and antenna pattern 
calculations. We will now present some preliminary results as a means of 
validating the numerical implementation. In particular, we shall compare 
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the numerically computed mode wavenumbers to the analytical solution for 
air-filled cavities. We shall also present some preliminary comparisons of 
RCS calculations with known results. Antenna pattern calculations will be 
deferred to a later report. 

7.1 Resonant Wavenumbers 

An important cavity characteristic is its resonant frequencies. These may be 
determined analytically from the eigenvalue series solution of the wave equa- 
tion for simple geometries. For example, if an air-filled rectangular cavity is 
considered, the resonant wavenumbers may be determined from the charac- 
teristic equation. Such a cavity is shown in figure 4a. Likewise, the resonant 
wavenumbers of a cylindrical sector cavity (see figure 4b) may be determined 
from its characteristic equation. The dimensions of the cylindrical sector cav- 
ity may be chosen so that this cavity is equivalent to a rectangular one which 
then requires that the resonant wavenumbers be identical. This was done for 
a equivalent rectangular cavity with dimensions: a = 0.25 cm, b = 0.5 cm, 
c = 1.0 cm. The following table illustrates the capability of the presented 
formulation to numerically compute the first five resonant wavenumbers 


mode 

exact 

computed 

error 

TE 1U 

5.236 

5.306 

1.34 

T M i io 

7.025 

7.441 

5.92 

T Equ 

7.531 

7.908 

5.01 

T Em 

7.531 

7.928 

5,27 

TMm 

8.179 

8.562 

4.68 


The mode number corresponds to the equivalent rectangular cavity and the 
same is true for the exact wavenumbers. The the computed wavenumbers 
were determined by solving the generalized eigenvalue problem (36) and we 
note that the percent error is comparable to the error reported using bricks 
to discretize a rectangular cavity [7]. 

A cylindrical sector cavity is shown in figure 4b. We have compared the 
exact resonant wavenumbers with computed results for a cavity where p b — 
5 cm, p a = 4.75 cm, a = <f> r - <j>i = 5°, z t = 0.5 cm, and z b = 0.0 cm in the 
following table of the first five modes. 
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mode 

exact 

computed 

error 

TE X u 

9.695 

9.856 

1.66 

T E\oi 

14.051 

14.283 

1.65 

TE n2 

14.575 

14.816 

1.65 

TM „o 

14.575 

15.539 

6.61 

TMiu 

15.872 

16.134 

1.65 


Again the percent error is within the expected bounds. The exact resonant 
wavenumbers of a cylindrical shell (see figure 4c) may also be found if the 
cavity is filled with a homogeneous material. For a shell where pi = 1.0 
cm, pi = 0.9 cm, z t = 1.0 cm and = 0.0 cm the comparison between the 

calculated and exact wavenumbers are given in the table below. 


mode 

exact 

computed 

error 

T M\oo 

6.246 

6.450 

3.27 

T Mi io 

6.393 

6.557 

2.57 

TE lu 

6.428 

6.608 

2.80 

T M\2o 

6.814 

6.713 

1.51 

T E\2i 

6.831 

6.733 

1.46 


7.2 Scattering 

The FEM-BI formulation has been shown to provide an excellent tool for 
computing the scattered field attributed to a cavity recessed in a ground 
plane [4]. This report presents an analogous method when the cavity is re- 
cessed in an infinite, metallic circular cylinder. We have shown that the finite 
element method may be used to numerically calculate interior resonances 
for closed cavities. We now wish to validate the proposed formulation for 
scattering calculations from an open cavity due to plane wave incidence. Un- 
fortunately, we currently do not have a suitable experimental or numerical 
data to determine the accuracy of this method. However, for small cavities 
embedded in a large radius cylinder, we expect our results to compare quite 
favorably with a similar method which was formulated for planar structures 
[4]. Below, we perform such comparisons for air-filled cavities. 

Comparisons have been made between a cavity recessed in a cylinder such 
as the one shown in figure 1 and an equivalent cavity in a metallic plane. First 
we consider the cavity whose dimensions are w x l x d — 2Ax2Ax 0.25A 


24 




cavity. The azimuthal width is given by w = where a is the cylinder 

radius in wavelengths and A(fi is the azimuthal angle subtended by the cavity. 
The axial length (/) and cavity depth (d) are normalized by the wavelength. 
The incident wave is given by (28) and we will consider only TM polarization 
where a = 0, and consequently E $ = 0. Figure 5 illustrates the comparison 
of the scattered fields for such a cavity in a cylinder with a = 20A and 
a similar cavity in a ground plane for normal incidence [0{ = 90°, <p x ~ 0°). 
Additional bistatic patterns are shown in figures 6 and 7 for oblique incidence 
with (6i = 45 °,<j>i = 0°) and (6 t = 20°, <f>i = 0°), respectively. In all of these 
three figures, the azimuth observation angle is <f> = 0° and the elevation angle 
(#) is allowed to vary. Figure 8 is a backscatter calculation of the same cavity 
for fixed <f> = 0° and varying 0. Finally, we repeat these calculations for a 
w x / x d = 2A x 0.5A x 0.25A cavity in figures 9-12. Clearly the two solutions 
are in excellent agreement. 

8 Future Work 

In this report, we have presented the formulation and preliminary validation 
of a FEM-BI method suitable for cavities recessed in a circular cylinder. 
Some validation of this formulation was presented for the special case of 
a large radius cylinder. However, we are actively seeking validation data 
where the curvature of the structure has a prominent effect on the scattering. 
Such validation data will either be found by measurement or comparison 
with another numerical solution. Such a solution will most likely be for an 
groove in a cylinder which will necessitate our consideration of a long narrow 
cavity and normal incidence only. Comparisons with curved cavities will help 
in determining which asymptotic Green’s function formula (18)-(20) is best 
suited for these calculations. It will also aid in the identification of the source 
causing the difficulties with the other principal polarization. Inaccuracies 
using Pathak’s formula appears to be the cause of this difficulty. 

After suitable scattering validation has been accomplished and a final 
decision as to which Green’s function formula is most accurate, we shall 
implement a host of features which are useful for antenna design. Metallic 
patches and resistive cards will be placed on the aperture and different feeds 
will be included. Since the mutual coupling between antenna array elements 
is expected to be crucial for conformal arrays, input impedance and mutual 
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impedance calculations will be also considered. 

Currently, we have implemented the formulation without exploiting the 
convolutional kernel of the boundary integral. In fact, we are using a sparse 
(direct) matrix solver requiring the full boundary integral matrix be stored in 
its entirety. Such a direct solver allows moderately large geometries without 
the need to resolve the system for each new incidence angle. However, for 
large problems, the boundary integral must be evaluated using a Fast Fourier 
Transform (FFT) in conjunction with an iterative solver such as the Bi- 
Conjugate Gradient (BiCG) method. The analogous FEM-BI formulation 
for planar structures [4] uses such a solver with impressive results both in 
terms of accuracy and efficiency. 
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Figure 10. Bistatic scattering by a 2A x 0.5A x 0.25A cavity in a cylinder 
for TM polarization and 0 t = 45 = <f) = 0° 
compared with an equivalent cavity recessed in a ground plane. 

Figure 11. Bistatic scattering by a 2 A x 0.5A x 0.25A cavity in a cylinder 
for TM polarization and 9, — 20 = <j> = 0° 
compared with an equivalent cavity recessed in a ground plane. 

Figure 12. Backscatter scattering by a 2A x 0.5A x 0.25A cavity in a cylinder 
for TM polarization and <f> — 0° compared with 
an equivalent cavity recessed in a ground plane. 
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Figure 1. Illustration of the cavity geometry situated 
on a metallic cylinder. 
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Figure 2. Cylindrical shell element. 






Figure 4. Cavity geometries for eigenvalue computation. 

a) Rectangular cavity. 

b) Circular cylinder sector cavity. 

c) Concentric circular cylinder cavity. 
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Figure 5. Bistatic scattering by a 2A x 2A x 0.25A cavity in a cylinder 
for Eg polarization and 0, = 90°, <f>{ = (j> — 0° 
compared with an equivalent cavity recessed in a ground plane. 
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Figure 6. Bistatic scattering by a 2A x 2A x 0.25A cavity in a cylinder 
for TM polarization and 0, = 45°, fa = <f> = 0° 
compared with an equivalent cavity recessed in a ground plane. 
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Figure 7. Bistatic scattering by a 2A x 2A x 0.25A cavity in a cylinder 
for TM polarization and 6 , = 20°, fa = <f> = 0° 
compared with an equivalent cavity recessed in a ground plane. 
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Figure 8. Backscatter scattering by a 2A x 2A x 0.25A cavity in a cylinder 
for TM polarization and <f> = 0° compared with 
an equivalent cavity recessed in a ground plane. 
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Figure 9. Bistatic scattering by a 2A x 0.5A x 0.25A cavity in a cylinder 
for TM polarization and 0, = 90°,^, = <f> = 0° 
compared with an equivalent cavity recessed in a ground plane. 
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Figure 10. Bistatic scattering by a 2A x 0.5A x 0.25A cavity in a cylinder 
for TM polarization and 0, = 45°, 4>i = (f> = 0° 
compared with an equivalent cavity recessed in a ground plane. 
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Figure 11. Bistatic scattering by a 2 A x 0.5A x 0.25A cavity in a cylinder 
for TM polarization and 0, = 20°, <f>i — (f> = 0° 
compared with an equivalent cavity recessed in a ground plane. 
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Figure 12. Backscatter scattering by a 2A x 0.5A x 0.25A cavity in a cylinder 
for TM polarization and <j> = 0° compared with 
an equivalent cavity recessed in a ground plane. 




